%=========================================================================%
% This program solves and estimates the parameters of the retirement 
% choice and asset accumulation model. There are three steps:
% 1. Obtain data moments (targets for estimation step)
% 2. Solve the model given a guess for the underlying parameters
% 3. Iterate step 2 by generating simulated versions of the data
%    moments and minimizing distance to the data moments. 
% 4. Obtain simulated retirement and saving histories in-simulation
%    using estimated parameters.
%=========================================================================%   

%Log on...
diary (sprintf('master_estimate_log_%s.txt',datestr(now,'yyyymmdd')));

%Code files (functions) are here
cd '/.../structural-model-code';
addpath ('...data/HRS MiCDA/Tables');
addpath ('...data/FSRDC/Tables');
addpath ('...data');
addpath ('.../structural-model-code/lininterp');

%-------------------------------------------------------------------------%
% Obtain data moments
%-------------------------------------------------------------------------%
% Treatment effect estimates for 56-64 y.o's starting from period 0 to 12
temp=csvread('te_emp_5664.csv',6,0,[6,0,18,1]);
TE_D=temp(1:end,1);
V_TE=temp(:,2).^2; %variances 
clear temp;

%LEHD control group employment rates
temp=csvread('Controllfp_levels.csv',1,0,[1,0,17,1]);
%Extract periods -5 to 12 
E_D=temp(:,1); 
V_E=temp(:,2).^2; %variances 
clear temp;

%Stack the moments (scale TE's)
scaler = 10;
dtarget{1,1} = [E_D;scaler.*TE_D];

%Use I as wt matrix  
dtarget{2,1} = eye(length(dtarget{1,1}));

%-------------------------------------------------------------------------%
% 2. Estimate parameters using SMM and  I as wt matrix 
% There are five parameters describing preferences
% i. sigma (EIS)
% ii. rho (AR(1) persistence of work disutility process)
% iii. sigma_v (AR(1) idiosyncratic variance of disutility process)
% iv. phi is the age trend slope
% v. gamma is the labor disutility constant
%-------------------------------------------------------------------------%

%Set seed
rng(1,'twister');

%-------------------------------------------------------------------------%
%Estimate parameters using first step weight matrix 
%-------------------------------------------------------------------------%
%Initialize parameters;
%sorted results of sobol sequence runs 
SMM_criteria = ...
    csvread('/shared/Patki_ra/Aging LFP/output/SMMcriteria_sobolseq_eis15_10wt.csv',0,0);
x0_sobol = round(SMM_criteria(1,2:5),4); 
f = @(p)smm_eis15(p,dtarget,scaler);
options = optimset('MaxIter',100,'Display','iter');
[x1,~,~,~] = fminsearch(f,x0_sobol,options)  
csvwrite('/shared/Patki_ra/Aging LFP/output/x1_eis15.csv',x1);

%obtain policy functions for estimated model
x1 = csvread('/shared/Patki_ra/Aging LFP/output/x1_eis15.csv');
[optfunc,optfunc_frz]=polfunc_function_eis15(x1);

%-------------------------------------------------------------------------%
% 3. Compute VCV of the simulated moments 
%-------------------------------------------------------------------------%
p_star1 = [1/(1+exp(x1(1)));exp(x1(2));x1(3);x1(4)];
%Define some input parameters:
rho = p_star1(1);
sigma_v = p_star1(2);
vcv_params=[rho;sigma_v];
%Omega is VCV
[Sm,Omega]=vcv_matrix(500,30,optfunc,optfunc_frz,vcv_params,scaler); 

%-------------------------------------------------------------------------%
% 4. Chi-2 model-fit statistics
%-------------------------------------------------------------------------%
%Data moments
DatM = dtarget{1,1};

%Moments evaluated at the parameter estimate (simulated moments)
[E_S,~,~,~,~,TE_S]=smom_function_eis15(x1); 
SimM = [E_S(:,1);TE_S];

theta1 = [1/(1+exp(x1(1)));exp(x1(2));x1(3);x1(4)];

%SMM objective function (J-stat)
m = SimM - DatM;
J = m'/Omega*m;
df = length(SimM)-length(theta1);
chi2test = [J df];

%-------------------------------------------------------------------------%
% 6. Output the results
%-------------------------------------------------------------------------%
csvwrite('.../model-output/chi2test_Iwt_eis15.csv',chi2test);
csvwrite('.../model-output/paramest_Iwt_eis15.csv',theta1);
%Simulated moments
csvwrite('.../model-output/E_S_eis15.csv',E_S);
csvwrite('.../model-output/TE_S_eis15.csv',TE_S);